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1  Introduction  and  motivation 

In  the  last  three  decades,  computer  simulation  tools  have  achieved  wide 
spread  use  in  the  design  and  analysis  of  engineering  devices.  This  has 
shortened  the  overall  product  design  cycle  and  it  has  also  provided  bet¬ 
ter  understanding  of  the  operating  behavior  of  the  systems  of  interest.  As 
a  consequence  numerical  simulations  have  lead  to  a  reduction  of  physical 
prototyping  and  to  lower  costs. 

In  spite  of  this  considerable  success,  it  remains  difficult  to  provide  objec¬ 
tive  confidence  levels  in  quantitative  information  obtained  from  numerical 
predictions.  The  complexity  arises  from  the  amount  of  uncertainties  related 
to  the  inputs  of  any  computation  attempting  to  represent  a  physical  system. 
As  a  result,  especially  in  the  area  of  reliability  and  safety,  physical  testing 
remains  the  dominant  mechanism  of  certification  of  new  devices.  Rigor¬ 
ous  quantification  of  the  errors  and  uncertainties1  introduced  in  numerical 
simulations  is  required  to  establish  objectively  their  predictive  capabilities. 

Procedures  to  establish  the  quality  of  numerical  simulations  have  been 
organized  within  the  framework  of  Verification  and  Validation  (V&V)  ac¬ 
tivities  [34],  Verification  is  a  mathematical  process  that  aims  at  answering 
the  question:  “are  we  solving  the  equations  correctly?”.  The  objective  is 
to  quantify  the  errors  associated  to  the  algorithms  employed  to  obtain  the 
solution  of  the  governing  equations.  Validation,  on  the  other  hand,  aims  at 
answering  the  question  “are  we  solving  the  correct  equations?”.  The  goal 
is  to  identify  the  appropriateness  of  the  selected  mathematical/physical  for¬ 
mulation  to  represent  the  device  to  be  analyzed.  Validation  always  involves 
comparisons  of  the  numerical  predictions  to  reality,  whereas  verification  only 
involves  numerical  analysis  and  tests. 

There  is  a  growing  recognition  of  the  fact  that  validation  cannot  be  car¬ 
ried  out  without  explicitly  accounting  for  the  uncertainties  present  in  both 
the  measurements  and  the  computations.  Experimentalists  are  typically  re¬ 
quired  to  report  uncertainty  bars  to  clearly  identify  the  repeatability  and  the 
errors  associated  to  the  measurements.  Validation  must  be  carried  out  ac¬ 
knowledging  the  nature  of  the  experimental  uncertainties  and  by  providing 
a  similar  indication  of  the  computational  error  bars.  One  of  the  objective 
of  uncertainty  quantification  methods  is  to  construct  a  framework  to  esti¬ 
mate  the  error  bars  associated  to  given  predictions.  Another  objective  is  to 
evaluate  the  likelihood  of  a  certain  outcome;  this  obviously  leads  to  better 
understanding  of  risks  and  improves  the  decision  making  process. 

1The  difference  between  errors  and  uncertainties  will  be  given  later 
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2  Definitions  and  basic  concepts 

The  uncertainty  quantification  community  has  introduced  precise  definitions 
to  characterize  various  types  of  uncertainties. 

2.1  Errors  vs.  uncertainties 

The  American  Institute  of  Aeronautics  and  Astronautics  (AIAA)  ’’Guide 
for  the  Verification  and  Validation  of  CFD  Simulations”  [35]  defines  errors 
as  recognisable  deficiencies  of  the  models  or  the  algorithms  employed  and 
uncertainties  as  a  potential  deficiency  that  is  due  to  lack  of  knowledge. 

This  definition  is  not  completely  satisfactory  because  does  not  precisely 
distinguish  between  the  mathematics  and  the  physics.  It  is  more  useful  to 
define  errors  as  associated  to  the  translation  of  a  mathematical  formulation 
into  a  numerical  algorithm  (and  a  computational  code). 

Errors  are  typically  also  further  classified  in  two  categories:  acknowl¬ 
edged  errors  are  known  to  be  present  but  their  effect  on  the  results  is  deemed 
negligible.  Examples  are  round-off  errors  and  limited  convergence  of  certain 
iterative  algorithms.  On  the  other  end,  unacknowledged  errors  are  not  rec¬ 
ognizable2  but  might  be  present;  implementation  mistakes  (bugs)  or  usage 
errors  can  only  be  characterized  by  comprehensive  verification  tests  and 
procedures. 

Using  the  present  definition  of  errors,  the  uncertainties  are  naturally  as¬ 
sociated  to  the  choice  of  the  physical  models  and  to  the  specification  of  the 
input  parameters  required  for  performing  the  analysis.  As  an  example,  nu¬ 
merical  simulations  require  the  precise  specification  of  boundary  conditions 
and  typically  only  limited  information  are  available  from  corresponding  ex¬ 
periments  and  observations.  Therefore  variability,  vagueness,  ambiguity  and 
confusion  are  all  factors  that  introduce  uncertainties  in  the  simulations.  A 
more  precise  characterization  is  based  on  the  distinction  in  aleatory  and 
epistemic  uncertainties. 

2.2  Aleatory  uncertainty 

Aleatory  uncertainty3  is  the  physical  variability  present  in  the  system  being 
analysed  or  its  environment.  It  is  not  strictly  due  to  a  lack  of  knowledge 

2In  principle,  using  the  AIAA  definition,  unacknowledged  errors  could  be  considered 
uncertainties  because  they  are  associated  to  lack  of  knowledge 

3Aleatory  uncertainty  is  also  referred  to  as  variability,  stochastic  uncertainty  or  irre¬ 
ducible  uncertainty. 
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and  cannot  be  reduced.  The  determination  of  material  properties  or  operat¬ 
ing  conditions  of  a  physical  system  typically  leads  to  aleatory  uncertainties; 
additional  experimental  characterization  might  provide  more  conclusive  ev¬ 
idence  of  the  variability  but  cannot  eliminate  it  completely.  Aleatory  uncer¬ 
tainty  is  normally  characterized  using  probabilistic  approaches. 

2.3  Epistemic  uncertainty 

Epistemic  uncertainty4,  is  what  is  indicated  in  the  AIAA  Guide  (AIAA 
1998)  as  ’’uncertainty”5,  i.e.  a  potential  deficiency  that  is  due  to  a  lack  of 
knowledge.  It  can  arise  from  assumptions  introduced  in  the  derivation  of 
the  mathematical  model  used  or  simplifications  related  to  the  correlation  or 
dependence  between  physical  processes.  It  is  obviously  possible  to  reduce  the 
epistemic  uncertainty  by  using,  for  example,  a  combination  of  calibration, 
inference  from  experimental  observations  and  improvement  of  the  physical 
models.  Epistemic  uncertainty  is  not  well  characterized  by  probabilistic  ap¬ 
proaches  because  it  might  be  difficult  to  infer  any  statistical  information 
due  to  the  nominal  lack  of  knowledge.  A  variety  of  approaches  have  been 
introduced  to  provide  a  more  suitable  framework  for  these  analysis.  Typi¬ 
cal  examples  of  sources  of  epistemic  uncertainties  are  turbulence  modeling 
assumptions  and  surrogate  chemical  kinetics  models. 

2.4  Sensitivity  vs.  uncertainty  analysis 

Sensitivity  analysis  (SA)  investigates  the  connection  between  inputs  and 
outputs  of  a  (computational)  model;  more  specifically,  it  allows  to  identify 
how  the  variability  in  an  output  quantity  of  interest  is  connected  to  an 
input  in  the  model  and  which  input  sources  will  dominate  the  response  of 
the  system.  On  the  other  hand,  uncertainty  analysis  aims  at  identifying  the 
overall  output  uncertainty  in  a  given  system.  The  main  difference  is  that 
sensitivity  analysis  does  not  require  input  data  uncertainty  characterization 
from  a  real  device;  it  can  be  conducted  purely  based  on  the  mathematical 
form  of  the  model.  As  a  conclusion  large  output  sensitivities  (identified 
using  SA)  do  not  necessarily  translate  in  important  uncertainties  because 
the  input  uncertainty  might  be  very  small  in  a  device  of  interest. 

SA  is  often  based  on  the  concept  of  sensitivity  derivatives,  the  gradient  of 
the  output  of  interest  with  respect  to  input  variables.  The  overall  sensitivity 

4Epistemic  uncertainty  is  also  called  reducible  uncertainty  or  incertitude 

5  Aleatory  uncertainty  is  not  mentioned  in  the  AIAA  Guidelines 
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is  then  evaluated  using  a  Taylor  series  expansion,  which,  to  first  order,  would 
be  equivalent  to  a  linear  relationship  between  inputs  and  outputs. 

3  Predictions  under  uncertainty 

Computer  simulations  of  an  engineering  device  are  performed  following  a 
sequence  of  steps. 

Initially  the  system  of  interest  and  desired  performance  measures  are 
defined.  The  geometrical  characterization  of  the  device,  its  operating  condi¬ 
tions,  the  physical  processes  involved  are  identified  and  their  relative  impor¬ 
tance  must  be  quantified.  It  is  worthwhile  to  point  out  that  the  definition  of 
the  system  response  of  interest  is  a  fundamental  aspect  of  this  phase.  The 
next  step  is  the  formulation  of  a  mathematical  representation  of  the  system. 
The  governing  equations  and  the  phenomenological  models  required  to  cap¬ 
ture  the  relevant  physical  processes  need  to  be  defined.  In  addition,  the  pre¬ 
cise  geometrical  definition  of  the  device  is  introduced.  This  step  introduces 
simplification  with  respect  to  the  real  system;  for  example  small  geometrical 
components  are  eliminated,  or  artificial  boundaries  are  introduced  to  reduce 
the  scope  of  the  analysis.  With  a  well  defined  mathematical  representation 
of  the  system,  the  next  step  if  to  formulate  a  discretized  representation. 
Numerical  methods  are  devised  to  convert  the  continuous  form  of  the  gov¬ 
erning  equations  into  an  algorithm  that  produces  the  solution.  This  step 
typically  requires,  for  example,  the  generation  computational  grid,  which  is 
a  tessellation  of  the  physical  domain.  Finally  the  analysis  can  be  carried 
out. 

The  introduction  of  uncertainty  in  numerical  simulations  does  not  alter 
this  process  but  introduces  considerable  difficulties  in  each  phase.  It  is  useful 
to  distinguish  three  steps:  data  assimilation,  uncertainty  propagation  and 
certification. 

3.1  Data  assimilation 

Data  assimilation  consists  of  a  study  of  the  system  of  interest  that  aims  at 
identifying  the  properties,  physical  processes  and  other  factors  required  to 
fully  characterize  it.  The  analysis  is  typically  focused  on  the  specific  in¬ 
puts  required  by  the  mathematical  framework  that  will  be  applied  in  the 
simulations.  As  an  example,  the  boundary  conditions  required  in  numeri¬ 
cal  simulations  should  be  inferred  from  observation  of  the  device  of  inter¬ 
est  or  specific  experiments.  Given  the  limited  degree  of  reproducibility  of 
experimental  measurements  and  the  errors  associated  to  the  measurement 
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techniques,  these  quantity  are  known  with  a  certain  degree  of  uncertainty 
(typically  specified  as  an  interval,  x  ±  y%).  Probabilistic  approaches  treat 
these  quantities,  that  overall  characterize  the  aleatory  uncertainty,  as  ran¬ 
dom  variables  assuming  values  within  specified  intervals.  In  mathematical 
terms  this  is  equivalent  to  define  random  variables  with  a  specified  probabil¬ 
ity  distribution  functions  (PDF)  [13].  The  obvious  choice  is  to  use  random 
variables  defined  using  analytical  distributions  (Gaussian,  uniform,  etc.).  It 
is  difficult  to  justify  this  choice  [5]  solely  from  experimental  evidence  be¬ 
cause  of  the  limited  amount  of  data  typically  available6;  in  many  situation 
the  only  data  available  are  obtained  from  expert  opinions  and  can  lead  to 
ambiguous  or  conflicting  estimates.  Alternative  approaches  have  been  de¬ 
vised  to  provide  a  more  flexible  framework  to  handle  this  situation,  evidence 
theory  is  one  such  approach  [37,  38]. 

In  the  context  of  probabilistic  approaches,  the  objective  of  data  assimila¬ 
tion  is  to  define  PDFs  of  each  of  the  input  quantity  used  in  the  computational 
tool. 

3.2  Uncertainty  propagation 

Once  probability  distributions  are  available  for  all  the  input  quantities  in 
the  computational  algorithm,  the  objective  is  to  compute  the  PDFs  of  the 
output  quantities  of  interest.  This  step  is  usually  the  most  complex  and 
computationally  intensive  for  realistic  engineering  simulations.  A  variety  of 
methods  are  available  in  the  literature,  from  sampling  based  approaches  (e.g. 
Monte  Carlo)  to  more  sophysticated  stochastic  spectral  Galerkin  approaches. 
In  the  next  section  these  methods  are  described  in  detail. 

3.3  Certification 

Once  the  statistics  of  the  quantity  of  interest  have  been  computed,  various 
metrics  can  be  used  to  characterize  the  system  output,  depending  on  the 
specific  application.  The  most  common  use  of  such  statistical  information 
is  a  reliability  assessments,  where  the  likelihood  of  a  certain  outcome  is  es¬ 
timated  and  compared  to  operating  margins.  In  a  validation  context,  the 
PDF  (or  more  typically  the  cumulative  distribution  function)  is  compared  to 
experimental  observation  to  extract  a  measure  of  the  confidence  in  the  nu¬ 
merical  results.  The  characterization  of  these  measures,  so-called  validation 

6Note  that  the  specification  of  an  interval  is  not  equivalent  to  a  specification  of  a 
uniform  probability  distribution.  An  interval  formally  does  not  contain  probabilistic  in¬ 
formation! 
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metrics,  is  an  active  area  of  research  [33] . 

4  Probabilistic  uncertainty  propagation 

Within  a  probabilistic  framework,  the  problem  of  uncertainty  propagation 
consists  of  the  generation  of  PDFs  of  the  outcomes  given  (known)  distri¬ 
bution  of  all  the  input  parameters.  Several  classes  of  methods  have  been 
developed  to  solve  this  problems;  in  this  section  three  popular  approaches 
are  described. 

Consider  the  vector  x  =  {x\ ,...xd)  containing  the  input  quantities  to 
the  computational  model;  assume  that  y  =  g(x)  is  the  output  of  interest;  g 
is  possibly  the  result  of  a  complex  fluid  dynamic  simulation. 

In  probabilistic  uncertainty  quantification  approaches  the  stochastic,  in¬ 
put  quantities  x  are  represented  as  independent  continuous  random  vari¬ 
ables  Xi(ui)7  mapping  the  sample  space  flj  to  real  numbers,  Xi  :  D*  — >  M. 
This  assumption  in  practical  terms  increases  the  dimensionality  of  the  prob¬ 
lem:  the  original  deterministic  outcome  y  =  f(xi,...,Xi,...xo)  becomes  a 
stochastic  quantity  y  =  f(xi,...,Xi,...xo  '■  uq, ...,  cjj, ...,  ud)-  The  objective 
is  to  compute  the  PDF  of  y,  fy,  in  order  to  evaluate  the  likelihood  of  a 
certain  outcome,  or,  in  general,  statistics  of  y.  The  expected  value  E[y]  and 
the  variance  Var[y\  are  defined  as: 


E[y]  =  /  zfy{z)dz  (1) 

J  OO 

POO 

V ar [y]  =  /  (z  -  E[y])2  fy(z)dz  =  E[y2]  -  (E[y])2  .  (2) 

J  OO 

Note  that  y  is  a  stochastic  variable  while  the  expected  value  and  the  variance 
are  deterministic  quantities. 

4.1  Sampling  techniques 

Sampling-based  techniques  are  the  simplest  approaches  to  propagate  un¬ 
certainty  in  numerical  simulations:  they  involve  repeated  simulations  (also 
called  realizations)  with  a  proper  selection  of  the  input  values.  All  the  results 
are  then  collected  to  generate  a  statistical  characterization  of  the  outcome. 
In  the  following  the  Monte  Carlo  approach  and  the  Latin  Hypercube  sam¬ 
pling  strategy  are  described. 

7The  assumption  of  independence  is  mathematically  stated  as  fxl,...,xD  =  fx1  ■  ■  •  fxD, 
where  fxl,...,xD  the  joint  PDF  is  the  product  of  the  marginals 
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Figure  1:  Equally  probability  intervals  defined  on  the  input  parameter  space 
for  a  Gaussian  PDF.  Each  of  the  interval  is  sampled  in  the  Latin  hypercube 
method. 


4.1.1  Monte  Carlo  method 

The  Monte  Carlo  method  [31]  is  the  oldest  and  most  popular  sampling  ap¬ 
proach.  It  involves  random  sampling  from  the  space  of  the  random  variables 
Xi  according  to  the  given  PDFs.  The  outcome  is  typically  organized  as  a 
histogram  and  the  statistics  are  readily  computed  from  Eq.  (2)  by  replac¬ 
ing  the  integrals  with  sums  over  the  number  of  samples.  The  method  has 
the  advantage  that  it  is  simple,  universally  applicable  and  does  not  require 
any  modification  to  the  available  (deterministic)  computational  tools.  It  is 
important  to  note  that  while  the  method  converges  to  the  exact  stochastic 
solution  as  the  number  of  samples  goes  to  infinity,  the  convergence  of  the 
mean  error  estimate  is  slow.  Hence  thousands  or  millions  of  data  samples 
may  be  required  to  obtained  accurate  estimations.  However,  the  conver¬ 
gence  does  not  directly  depend  on  the  number  of  random  variables  in  the 
problem.  In  this  form  the  Monte  Carlo  methods  always  give  the  correct 
answer,  but  a  prohibitively  large  number  of  realizations  may  be  required  to 
accurately  estimate  responses  that  have  a  small  probability  of  occurrence. 
On  the  other  hand,  the  convergence  of  the  low  order  statistics  (expected 
value  and  variance)  require  much  smaller  number  of  samples. 
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Figure  2:  Sampling  based  techniques  with  two  uncertain  inputs,  (left)  Monte 
Carlo;  (center)  Latin  Hypercube;  (right)  Lattice  based 


4.1.2  Latin  Hypecube  approach 

Several  methods  have  been  developed  to  accelerate  the  Monte  Carlo  ap¬ 
proach.  One  of  the  most  successful  is  the  Latin  Hypercube  sampling  (LHS) 
[29]  approach.  The  range  of  each  input  random  variable  X{  is  divided  in 
M  intervals  with  equal  probability.  In  Fig.  1  this  is  illustrated  for  a  single 
input  Gaussian  random  variable.  M  random  samples  are  drawn,  one  from 
each  of  the  intervals  and  the  convergence  is  faster  than  Monte  Carlo  because 
the  occurrence  of  low  probability  samples  is  reduced.  The  construction  is 
especially  effective  in  multiple  dimensions;  in  Fig.  2  an  example  of  samples 
generated  using  Monte  Carlo  and  Latin  Hypercube  for  two  uniform  ran¬ 
dom  variables  is  reported.  The  LHS  samples  projected  on  each  of  the  axis 
provide  optimal  coverage  of  the  parameter  space  [30].  On  the  other  hand, 
the  Monte  Carlo  samples  show  clusters  and  holes  that  eventually  lead  to 
the  slower  convergence.  In  Fig.  2  another  strategy  is  illustrated:  Lattice- 
based  sampling.  In  this  approach,  the  parameter  space  is  discretized  using 
regularly  spaced  points  (lattice);  once  the  solution  is  evaluated  at  those  lo¬ 
cations  a  random  shift  is  applied  to  the  lattice  and  another  set  of  solutions 
are  computed[31].  The  choice  of  the  lattice  is  related  to  the  distribution  of 
the  input  variables  similarly  to  the  LHS  approach  and  the  objective  is  again 
to  reduce  the  occurrence  of  clusters  and  holes  in  the  samples. 

4.2  Quadrature  methods 

One  of  the  objective  of  UQ  propagation  methods  is  the  computation  of 
the  statistics  of  an  outcome  of  interest,  such  as  its  expectation  and  the 
variance.  As  shown  earlier,  these  require  the  evaluation  of  integrals  (over 
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the  parameter  space)  and  it  is  therefore  natural  to  employ  conventional 
numerical  integration  techniques* 8. 

Let’s  consider  a  problem  with  one  uncertain  parameter  £;  the  objective 
is  to  compute  integrals  of  y(£).  A  class  of  quadrature  rules  are  based  on  in¬ 
terpolating  basis  functions  that  are  easy  to  integrate,  typically  polynomials. 
The  integral  is  expressed  as  a  weighted  sum  of  the  integrand  y  evaluated 
in  a  finite  number  of  locations  on  the  £-axis  (abscissas).  The  choice  of  the 
polynomial  basis  defines  the  weights  and  the  corresponding  abscissas. 

The  simplest  example  is  the  midpoint  rule  in  which  a  sequence  of  M 
equally  spaced  values  in  the  interval  [a,  b ]  (x\  =  a  and  xm  =  b)  is 
considered.  The  function  y  is  evaluated  at  M  —  1  abscissas  as  y((£W  -f 
£(*+!)) /2)  =  yi+l/2 .  The  integral  is  simply  obtained  as  the  sum  of  the  M  —  1 
values  |_i /2  times  the  size  of  the  interval  (b  —  a).  In  this  case  the  basis 
functions  are  piecewise  constant  and  the  weights  are  all  equal  to  ( b  —  a)/M. 
Quadratures  based  on  equally  spaced  abscissas  include  the  commonly  used 
trapezoidal  and  Simpson  rules  and  are,  in  general,  referred  to  as  Newton- 
Cotes  formulas. 

4.2.1  Stochastic  collocation 

Stochastic  collocation  refers  to  quadrature  methods  used  to  compute  in¬ 
tegrands  of  random  variables,  thus  over  a  stochastic  domain.  Although 
Newton-Cotes  formulas  are  applicable  in  this  context,  it  is  usually  prefer¬ 
able  to  consider  more  general  approaches,  in  which  the  abscissas  are  not 
equally  spaced.  Gaussian  quadratures  are  popular  in  the  field  of  uncertainty 
analysis  because  of  their  high  accuracy9;  the  general  format  is: 

rb  M 

/  y(t,)f(Od£  =  ^2w(tk)y(£k)  +  Rm{v )  (3) 

k= 0 

where  the  abscissas  £*.  are  associated  with  zeros  of  orthogonal  polynomials 
and  Wk  are  weights  and  /(£)  is  a  weighting  function;  in  the  present  context, 
/  represents  a  PDF.  Rm  is  the  remainder  term  which  determines  the  order 
of  accuracy  of  the  integration.  In  particular,  the  formula  (3)  is  said  to  have 
polynomial  degree  of  exactness  p  if  R.m(jj)  =  0  when  y(£)  is  a  polynomial 
of  degree  <  p.  Newton-Cotes  formulae  have  p  =  M  —  1  whereas  Gaussian 

sNumerical  integration  and  quadrature  are  typically  synonym  for  one  dimensional  func¬ 

tions.  In  higher  dimensions  the  term  cubature  is  often  used  instead  of  quadrature 

9  High  order  polynomial  interpolating  functions  constructed  on  equally  spaced  grids 
also  suffer  from  the  Runge’s  phenomenon. 
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Figure  3:  (Legendre  (left)  and  Hermite  (right)  polynomials. 


quadratures  achieve  higher  p  with  the  same  number  of  function  evaluations 
(M  +  1),  specifically  p  =  2 M  +  1.  This  is  the  main  advantage  of  Guassian 
quadrature  and  is  accomplished  by  virtue  of  the  orthogonality  condition  of 
the  interpolating  basis  functions [24]  [4],  which  is  in  general  stated  as: 

{Pm,Pn)  =  [  Pm(OPn{Of(Od£  =  Snm  (4) 

J  a 

where  5nm  is  the  Kroneker  symbol. 

The  most  commonly  used  form  of  Guassian  quadrature  is  the  Gauss- 
Legendre  integration  formula  which  is  based  on  Legendre  polynomials  4  (0  ■ 
In  Fig.  3  the  first  few  polynomials  are  reported.  These  are  constructed 
using  the  three  term  recurrence  relation  (with  £q(£)  =  1  and  4(£)  =  £): 


(n  +  1)4+1  (£)  =  (2ri  +  1)£4(0  -  tc4-i(£) 

The  weight  function  in  Eq.  (3)  is  simply  /(£)  =  1;  the  integration 
weights  are  computed  as 

Wk  =  f 

J  a 

and  the  abscissas  4  for  an  M-point  quadrature  rule  are  the  zeros  of  the 
£m+ i  polynomial.  As  an  example  the  computation  of  expectation  of  an 
output  in  the  presence  of  an  uncertain  input  defined  in  a  interval  [a,  b\  (as 
a  uniform  random  variable)  can  be  computed  as: 

rb  M 

E(y)  =  /  y(€)fy(€)d£  =  ^2y(£k)wk  (5) 

k= 0 
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For  uncertain  parameters  defined  on  unbounded  domains  such  as  ran¬ 
dom  variables  described  by  Gaussian  probability  distributions,  the  Gauss- 
Hermite  quadrature  is  used  [19].  In  this  case  the  three  term  recurrence 
relation  is  (/io(£)  =  1  and  /ii(£)  =  £): 

K+ 1(0  =  2£/in,(£)  -  2 nhn-x(£) 

and  the  weight  function  is  /(£)  =  e~^2 .  In  this  case  the  integration  weights 
Wk  are  expressed  as: 


/OO 

hk(£)e-?dZ 

-oo 

and  the  expectation  can  be  computed  as  before  just  evaluating  the  outcome 
at  M  +  1  abscissas  that  correspond  to  the  zeros  of  the  Hermite  polynomial 
h>M+ 1- 

There  is  an  obvious  connection  between  the  choice  of  the  polynomial 
and  the  probability  distribution  of  the  input  variables:  the  orthogonality 
condition  (Eq.  4).  In  the  present  context,  in  order  to  achieve  the  polyno¬ 
mial  exactness  of  Gaussian  quadratures  for  computing  the  expectation,  the 
polynomials  have  to  be  orthogonal  with  respect  to  the  weight  function  that 
corresponds  to  the  PDF  of  the  input  variable. 

The  two  Gauss  quadrature  rules  described  above  require  the  computation 
of  a  new  set  of  abscissas  (and  corresponding  function  evaluations)  for  each 
desired  order  M ;  alternative  approaches  have  the  property  that  the  abscissas 
are  nested,  e.g.  different  accuracy  formulas  share  abscissas.  A  popular 
nested  formulation  rule  is  the  Clenshaw-Curtis  quadrature[17]. 

In  practical  terms,  collocation  methods  for  uncertainty  propagation  re¬ 
quire  the  evaluation  of  zeros  and  weights  for  a  family  of  orthogonal  basis 
functions;  these  can  be  computed  and  stored  in  advance.  A  set  of  M  +  1 
independent  computations  are  performed  and  the  results  are  combined  to 
obtain  the  statistics  of  the  output  of  interest.  Collocation  can  therefore 
be  interpreted  as  a  sampling  technique;  it  retains  the  main  advantage  of 
the  Monte  Carlo  method  because  it  does  not  require  modifications  to  the 
existing  computational  tool. 

The  evaluation  of  the  PDF  of  the  output  quantities  is  somewhat  more 
complicated  for  stochastic  collocation  methods  than  the  computations  of 
the  output  statistics.  The  first  step  is  the  construction  of  the  polynomial 
interpolant  of  the  solution  y(£)  in  the  parameter  space  which  involves  the 
use  of  the  M  + 1  solutions  y(£W).  At  this  point,  the  interpolant  is  used  as  a 
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Figure  4:  Tensor  (left)  and  sparse  (right)  grid  quadrature  rules  for  two- 
dimensional  functions. 

replacement  for  the  original  function10  and  Monte  Carlo  sampling  is  used. 

4.2.2  Extension  to  multiple  dimensions 

The  natural  extension  of  interpolation  techniques  to  multiple  dimensions 
is  a  tensor  product  of  one-dimensional  interpolants.  Let  G(£i,  •  •  • ,  £d)  be 
the  outcome  of  interest,  expressed  as  a  function  of  D  random  variables. 
The  set  of  points  for  the  multi-dimensional  interpolant  is  the  tensor  prod¬ 
uct  of  the  abscissas  for  the  one-dimensional  interpolants;  in  the  case  that 
the  polynomial  order  is  the  same  in  every  direction,  the  total  number  of 
abscissas  is  (M  +  1  )D.  Although  the  tensor  product  extension  of  the  one- 
dimensional  quadrature  rules  is  simple  and  accurate,  it  clearly  suffers  from 
the  exponential  increase  in  the  number  of  abscissas  -  and  consequently  re¬ 
quired  function  evaluations  -  as  the  number  of  dimensions  increases  [20] .  To 
combat  this,  some  have  proposed  using  an  alternative  extension  based  on 
the  sparse  grid  construction,  attributed  to  Smolyak  [15].  In  this  case  the 
multi-dimensional  product  is  constructed  as  a  linear  combination  of  tensor 
product  interpolants,  each  with  relatively  few  points  in  its  respective  point 
set.  This  greatly  reduces  the  total  number  of  abscissas  -  particularly  as 
the  dimension  D  increases  -  while  retaining  the  accuracy  and  convergence 

10  The  approximated  function  is  typically  called  a  surrogate  response  surface 
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properties  of  one-dimensional  interpolants. 

Fig.  4  illustrates  the  extension  of  Gauss-Legendre  quadrature  in  two- 
dimensions  using  tensor  products  and  sparse  grid  rules.  Different  set  of 
abscissas  (M  is  varied  from  1  to  11)  are  reported  to  illustrate  how,  for 
sparse  grids,  the  number  of  points  and  corresponding  function  evaluations 
remains  relatively  low  even  for  high  order  polynomial  interpolants. 

Sparse  grid  constructions  reduce  the  curse  of  dimensionality  and  con¬ 
verge  very  rapidly  if  the  function  of  interest  is  smooth.  Precise  convergence 
results  and  comparisons  to  Monte  Carlo  methods  are  reported  in  [9].  Exten¬ 
sions  of  this  approach  to  construct  anisotropic  grids  are  currently  an  active 
area  of  research;  the  objective  in  this  case  is  to  increase  the  order  of  the 
interpolating  polynomials  only  in  selected  directions  to  further  reduce  the 
computational  cost. 

Other  integration  methods  in  high  dimensions  are  also  available  in  the 
literature,  for  example  Stroud  formulae  [10]  although  have  not  yet  been 
explored  in  the  area  of  uncertainty  quantification. 

4.3  Spectral  methods 

The  third  class  of  methods  for  uncertainty  propagation  is  based  on  spectral 
methods  in  which  the  solution  is  expressed  using  a  suitable  series  expansion. 
These  approaches  are  intrusive  in  the  sense  that  the  mathematical  formula¬ 
tion  requires  modification  of  the  existing  deterministic  codes  -  in  contrast  to 
the  sampling  approaches  and  the  stochastic  collocation  methods  described 
before.  For  this  reason,  spectral  methods  have  to  be  described  with  reference 
to  a  given  mathematical  problem;  in  the  following  the  Burgers  equation  is 
used. 

4.3.1  Stochastic  Galerkin  approach 

In  this  approach  the  unknown  quantities  are  expresses  as  an  infinite  series  of 
orthogonal  polynomials  defined  on  the  space  of  he  random  input  variable11. 
This  representation  has  its  roots  in  the  work  of  Wiener  [23]  who  expressed  a 
Gaussian  process  as  an  infinite  series  of  Hermite  polynomials.  Ghanem  and 
Spanos  [22]  truncated  Wiener’s  representation  and  used  the  resulting  finite 
series  as  a  key  ingredient  in  a  stochastic  finite  element  method. 

As  an  example,  consider  the  ID  Burgers  equations 

ut  +  uux  =  0 

11Stochastic  Galerkin  methods  based  on  polynomial  expansions  are  often  referred  to  as 
Polynomial  Choas  approaches  [22] 
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and  assume  that  the  initial  condition  is  uncertain  and  characterized  by  one 
parameter,  £  a  Guassian  random  variable. 

u(x,  t  =  0)  =  £  *  sin(x) 

In  stochastic  Galerkin  approaches  the  (unknown)  solution  u  =  u(x,t ,  £) 
is  expressed  as  a  spectral  expansion  as 

OO 

u(x,t,Q  =  ^Tui(x,t)'I?i(€)  (6) 

i= 0 

where  the  coefficients  Ui(x,t )  are  independent  of  the  random  variable  and, 
therefore,  deterministic  (they  are  still  function  of  the  physical  dimensions). 
The  4q(£)  are  Hermite  polynomials  and  form  a  complete  orthogonal  set  of 
basis  functions. 

The  expansion  (6)  is  truncated  to  a  certain  order  M  and  the  approxi¬ 
mated  expression  for  u  is  inserted  in  the  governing  equation: 


Multiplying  by  T^(^)  and  integrating  over  the  probability  space  -  Galerkin 
projection  -  we  obtain  a  system  of  M  coupled  &  deterministic  equations  for 
the  coefficients  uf 


Mr,  M  M  r. 

E  +  =  0 


i= 0 


i= 0  j= 0 


for  fc  =  0, 1, ...,  M.  (7) 


Note  that  in  Eq.  (7)  the  double  and  triple  products  of  the  Hermite 
polynomials  are  simply  numbers  [24]  (cfr.  Eq.  4): 


(*iVj)  =  Sijil 


('kjTj'kfc)  =0  if  i  +  j  +  k  is  odd  or  max(i,  j.  k)  >  s 


(ViVjVk) 


i\j\k\ 

(s  —  i)\(s  —  j)!(s  —  k)! 


otherwise 


and  s  =  (i  +  j  +  k)/ 2. 

The  statistics  of  the  solutions  can  be  readily  computed  from  the  expansions 
coefficients.  In  particular  the  expectation  is 
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E[m]  =  E 


M 


y 


.k- 0 


M 


=  mqE^o]  +  ^2  UkK[uk}  =  U0 


fc= 1 


for  the  orthogonality  of  the  polynomial  basis.  The  variance  is  computed  as: 


Var[u]  =  E  [(m  -  E[m])2]  =  E 


E 


M  M 

X>*e  [*2]  =E^(fc!)2 


fc= i 


fc= i 


As  an  example,  the  polynomial  chaos  formulations  for  the  Burgers  equations 
can  be  written  explicitly.  For  M  =  2  the  expansion  (6)  is  m(x,  t,  £)  =  up  +  £mi  and 
the  governing  system  of  equations  is: 


du0 

dt 

du\ 

~dt 


Up 


+  Ul 


dup 

dx 

dup 

dx 


■  u  1 


Up- 


dui 

dx 

du\ 

dx 


0 

0 


For  M  =  3  the  expansion  (6)  is  u(x,t,£)  =  up  +  +  2(^2  —  1  )u2  and  the 

governing  system  of  equations  is: 


du  i 
~dt  "* 
du2 
~dt 


u  1 


dwo 

dt 

dup  _ 
dx 
dup 


du  o  chzi 
- 1" 

dx 
du  i 
dx 


U2 


dx 
{up  +  2m2) 
9m  i 


9a; 


■  Mr 


9a; 


2m2 


2mi 


(m0  +  4m2) 


9m2 

9a; 

9m2 

9a; 

9m2 

9a; 


0 

0 

0 


It  is  worth  to  note  that  the  choice  of  the  polynomial  basis  is  again  connected 
to  the  distribution  of  the  input  random  variables.  As  for  stochastic  collocation 
approaches  Hermite  and  Legendre  polynomials  are  used  for  Gaussian  and  uniform 
random  variables  respectively.  A  general  formulation  based  on  the  Askey  polyno¬ 
mials  is  presented  in  [21]. 

While  the  stochastic  Galerkin  technique  has  been  shown  to  produce  highly  ac¬ 
curate  statistics,  it  does  requires  extensive  modifications  to  existing  computational 
methods  to  solve  for  the  coefficients  of  the  expansion. 


4.3.2  Extension  to  multiple  dimensions 

The  presence  of  multiple  uncertain  parameters  increases  considerably  the  computa¬ 
tional  complexity  of  the  stochastic  Galerkin  approach  although  it  does  not  change 
the  overall  formulation. 
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It  is  formally  possible  to  formulate  the  basis  polynomials  in  terms  of  a  full 
tensor  product  of  one-dimensional  polynomials  similarly  to  what  considered  for  the 
stochastic  collocation  approach.  Ghanem  and  Spanos  [22]  introduced  a  different 
expansion.  For  example  for  a  problem  with  two  uncertain  parameters  expressed  in 
terms  of  Gaussian  random  variables  the  first  few  bi-variate  Hermite  polynomials 
are: 


•hotel  £2)  =  1 

•hitete)  =  Ci 
•t^tete)  =  £2 
•hstete)  =  £1-1 
$4tete)  =  C1C2 
^stete)  =  Cl  -  1 


(8) 


In  multiple  dimensions  the  number  of  coefficients  M  is  determined  by  the  num¬ 
ber  of  dimensions  D  and  the  highest  polynomial  order  P : 


M  = 


(. D  +  P)\ 
D\P\ 


(9) 


As  in  stochastic  collocation  approach,  the  stochastic  Galerkin  method  suffers 
from  the  curse  of  dimensionality.  Smolyak-type  constructions  have  not  been  at¬ 
tempted  for  the  polynomial  basis  used  in  the  spectral  expansions  but  offer  a  possible 
remedy. 


5  Examples 

The  examples  reported  in  the  following  offer  indications  of  the  relative  computa¬ 
tional  cost  of  the  uncertainty  propagation  approaches  illustrated  before.  In  these 
test  cases  the  uncertainty  is  assumed  as  part  of  the  problem  specification. 


5.1  Steady  Burgers  equation 

The  first  test  problem  is  based  on  the  ID  viscous  Burgers  equations  [36]: 


1  du  d2u 

2  (1  “  dx  = 

where  the  convective  flux  is  modified  so  that  an  exact  solution  can  be  obtained 
(manufactured  solution  [34]): 


u(x) 


1 

2 


1  +  tank 
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Figure  5:  Expectation  of  the  solution  for  the  Burgers  steady  state  problem 
with  uncertain  viscosity.  Monte  Carlo  sampling  with  different  number  of 
realization.  A  numerical  solution  computed  using  32  cells  (left)  is  compared 
to  the  exact  solution  (right). 


In  this  case  it  is  assumed  that  the  viscosity  is  uncertain,  for  example  as  a 
consequence  to  an  unknown  fluid  temperature.  A  Gaussian  random  variable  is 
considered  with  mean  and  variance  given  by  E[/r]  =  0.25  and  var[g s]  =  0.002512 
Calculations  have  been  performed  using  a  finite  volume  discretization  on  the 
domain  [-3:3];  boundary  conditions  are  derived  from  the  exact  solution.  Monte 
Carlo  sampling  is  initially  used  as  the  uncertainty  propagation  technique.  The 
expectation  of  the  velocity  is  reported  in  Fig.  5.1  and  has  been  computed  using 
a  different  number  of  samples  ranging  from  10  to  10,000.  The  results  shows  only 
small  dependency  on  the  number  of  realization;  in  the  same  figure  the  statistical 
convergence  is  also  evaluated  using  the  exact  solution  of  the  deterministic  problem. 
In  Fig.  5.1  the  solution  variance  is  reported.  In  this  case  it  is  clear  that  10  samples 
are  not  sufficient  and  converged  results  are  only  obtained  for  1,000  and  10,000 
samples.  Note  that  the  variance  is  zero  at  x  =  0  to  satisfy  the  symmetry  constraint. 

As  a  second  step  the  stochastic  Galerkin  method  has  been  implemented  and 
tested.  In  this  case  the  expansion  is  based  on  Hermite  polynomials.  Computations 
have  been  performed  using  the  same  discretization  used  in  the  deterministic  code 
[36].  Fig.  5.1  shows  results  obtained  using  different  order  of  expansions,  from  1  to  3; 
the  comparisons  to  Monte  Carlo  show  that  for  this  simple  problem  -  characterized  by 
small  variability  -  even  the  low  order  expansion  provide  a  reasonable  representation 

12It  should  be  noted  that  a  Gaussian  distribution  for  the  viscosity  implies  that  there  is 
a  non-zero  probability  of  having  very  low  -  or  even  negative  -  values  of  the  viscosity.  In 
this  case  the  high  mean  value  and  the  small  variance  reduce  the  occurrence  of  negative 
viscosity  to  virtually  zero. 
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Figure  6:  Variance  of  the  solution  for  the  Burgers  steady  state  problem 
with  uncertain  viscosity.  Monte  Carlo  sampling  with  different  number  of 
realization.  A  numerical  solution  computed  using  32  cells  (left)  is  compared 
to  the  exct  solution  (right). 
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Figure  7:  Expectation  (left)  and  variance  (right)  of  the  solution  for  the  Burg¬ 
ers  steady  state  problem  with  uncertain  viscosity.  The  Stochastic  Galerkin 
method  is  used  and  the  results  are  compared  to  Monte  Carlo  sampling. 
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Figure  8:  Computational  grid  for  the  flow  around  an  array  of  cylinders  with 
uncertain  surface  heat  flux 


of  the  variance  of  the  solution. 

5.2  Heat  transfer  computations  under  uncertainty 

In  this  problem  the  goal  is  the  computation  of  the  wall  temperature  for  an  array  of 
cylinders  in  crossflow  under  uncertain  free  stream  conditions  and  surface  heat  flux. 
The  mathematical  framework  is  based  on  the  two-dimensional  Reynolds-averaged 
Navier-Stokes  equations  and  the  computational  domain  is  a  channel  with  a  specified 
inflow,  periodic  boundary  conditions  on  the  top  and  bottom  boundaries,  and  an 
outflow.  An  unstructured  grid  is  used,  see  Fig.  5.2. 

We  assume  that  the  sources  of  uncertainty  are  the  specification  of  the  velocity 
boundary  condition  on  the  incoming  flow  and  the  definition  of  the  thermal  condition 
on  the  surface  of  the  cylinder.  The  inlet  velocity  profile  is  constructed  as  a  linear 
combination  of  two  cosine  functions: 

Uiniet{y)  =  1  +  0.25(6  cos(2t ry)  +  6  cos(10t ry)).  (10) 

where  6  and  6  are  two  uniform  random  variables  in  the  interval  [-1:1];  this  inflow 
specification  ensures  that  the  amplitude  of  the  random  fluctuations  did  not  cause 
the  velocity  to  be  negative.  The  wave  numbers  2  and  10  were  chosen  to  give  one 
low  and  one  high  frequency  fluctuation  while  maintaining  symmetry  in  the  problem. 
The  heat  flux  is  specified  as  an  exponential  function  of  6  over  the  cylinder,  namely 

(X)  =  e  — (°.l+0.05C3)(z— 0.5)  (11) 

OTl  cylinder 

where  n  is  the  normal  to  the  cylinder  surface.  The  heat  flux  is  slightly  greater 
at  the  left  side  of  the  cylinder  at  the  stagnation  point;  the  value  of  £3  determines 
precisely  how  much  greater.  Again  6  is  a  uniform  random  variable  in  [-1:1]. 

We  employed  Monte  Carlo  sampling  and  stochastic  collocation  to  solve  this 
problem.  In  particular,  we  used  both  tensor  product  and  sparse  grid  formalism  to 
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l 

#  tensor  grid  points 

#  sparse  grid  points 

2 

27 

25 

3 

125 

111 

4 

729 

351 

5 

2187 

1026 

Table  1:  Number  of  function  evaluation  required  to  achieve  a  certain  degree 
of  polynomial  exactness  using  tensor  product  and  sparse  grid  formulations. 
The  parameter  l  is  called  the  level  and  is  directly  connected  to  the  accuracy 
of  the  interpolant  on  tensor  grids;  on  sparse  grids  the  relationship  is  more 
complex. 


construct  quadrature  rules  in  the  three-dimensional  parameter  space  (£1,  £2,  £3).  In 
Table  1  the  number  of  function  evaluations  (abscissas)  required  to  achieve  a  certain 
order  of  polynomial  exactness  are  reported. 

The  expectation  and  variance  of  the  wall  temperature  are  reported  in  Fig.  5.2 
and  5.2,  respectively. 

In  this  application,  it  is  clear  that  accurate  temperature  statistics  can  be  ob¬ 
tained  efficiently  using  stochastic  collocation  with  respect  to  Monte  Carlo  sampling. 
The  difference  between  tensor  product  and  sparse  grid  formulations  is  very  small. 


5.3  Shock  dominated  problems 

Uncertainty  propagation  algorithms  based  on  polynomial  basis  are  particularly  ef¬ 
fective  when  the  output  of  interest  are  smoothly  varying  with  respect  to  the  input 
uncertainties.  The  presence  of  discontinuities  in  the  response  surface  poses  a  chal¬ 
lenge;  in  physical  terms,  these  discontinuities  represent  sharp  system  transitions 
such  as  the  occurrence  of  flow  separation  or  shock  interactions. 

Consider  the  problem  governed  by  the  non- homogenous  Burgers  equations  [27] : 


du  d  u2 
dt  dx  2 


d  sin2  x 
dx  2 


0  <  £  <  7T,  t  >  0 


with  the  initial  condition  u(x,  0)  =  /3sin£,  and  boundary  conditions  u(0,t)  = 
u(Tr,t)  =  0.  (3  represents  the  uncertainty  in  the  initial  condition  and  is  defined 
as  _ 


-1  +  v/l  +  4<r2£ 
2  ct2£ 


where  £  is  a  Gaussian  random  variable  and  a  is  used  to  control  the  amount  of 
uncertainty. 

The  exact  steady  solution  contains  a  shock  located  at  xs  =  /(/?)  and  the  desired 
output  is  the  PDF  of  the  shock  location.  The  response  surface  u  as  a  function  of  [3 
is  reported  in  Fig.  12. 
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Figure  9:  Expectation  of  the  wall  temperature.  Tensor  product  (left)  ver¬ 
sus  sparse  grids  (right).  Monte  Carlo  (MC)  statistics  are  based  on  10,000 
samples;  multiple  quadrature  formulas  are  reported  for  different  levels  (see 
Table  1). 


Figure  10:  Variance  of  the  wall  temperature.  Tensor  product  (left)  versus 
sparse  grids  (right).  Monte  Carlo  (MC)  statistics  are  based  on  10,000  sam¬ 
ples;  multiple  quadrature  formulas  are  reported  for  different  levels  (see  Table 
1)- 
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Figure  11:  Probability  distribution  function  of  the  wall  temperature  at  the 
stagnation  point.  Tensor  product  (left)  versus  sparse  grids  (right).  Monte 
Carlo  (MC)  statistics  are  based  on  10,000  samples;  multiple  quadrature 
formulas  are  reported  for  different  levels  (see  Table  1). 


Figure  12:  Solution  of  the  in-homogeneous  Burgers  equations  with  uncertain 
initial  condition.  The  sharp  transition  corresponds  to  a  shock  whose  location 
is  a  function  of  the  input  uncertainty. 
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Figure  13:  Probability  distribution  of  the  shock  location  for  the  inhomo¬ 
geneous  Burgers  equations  with  uncertain  initial  condition.  Low  ,  a  =  0.1 
(left)  and  high  input  variability,  a  =  0.6  (right).  Results  obtained  using 
high-order  stochastic  Galerkin  expansions  [28] 


Hestevan  et  al.  [25]  have  applied  the  stochastic  Galerkin  approach  to  study 
this  problem;  high  order  expansions  in  Hermite  polynomials  were  used  to  compute 
the  PDF  of  the  shock  location.  The  availability  of  exact  PDFs  allows  to  identify 
the  limitation  of  the  numerical  approach.  In  Fig.  5.3  the  solution  is  reported  for 
low  and  high  degree  of  input  variability  ( a  set  to  0.1  and  0.6  respectively).  For 
the  low  variability  case  the  computed  results  are  in  reasonable  agreement  with 
the  exact  PDF,  but  some  oscillations  are  present.  For  the  high  variability  case, 
the  oscillations  are  the  dominant  feature  and  the  shape  of  the  PDF  is  completely 
misrepresented.  In  this  case,  the  problem  is  that  the  solution  is  discontinuous  and 
the  polynomial  basis  (underlying  the  stochastic  expansion)  is  highly  oscillatory  - 
Gibbs  phenomenon. 

An  alternative  response  reconstruction  has  been  developed  to  handle  strong 
non-linear  or  discontinuous  surfaces.  The  method  is  based  on  a  Pade-Legendre 
interpolatory  formula  [32];  it  is  a  stochastic  collocation  approach  that  eliminates 
the  Gibbs  phenomenon.  Assume  N  solution  evaluations  are  available;  given  the 
integers  M  and  L,  the  pair  of  Legendre  polynomials  P  and  Q  of  order  less  or 
equal  than  M  and  L,  respectively,  are  said  to  be  the  solution  of  the  (N,  M,  L ) 
Pade-Legendre  interpolation  problem  [26]  of  u  if 

(P-Qu,^)=0  (12) 

and  the  rational  function  R(u)  :=  P/Q  is  defined  as  the  approximation  of  u.  The 
details  of  the  construction  of  R  are  reported  in  Chantrasmi  et  al.  [32];  here  it  is 
sufficient  to  observe  that  if  u  is  discontinuous,  the  polynomial  Q  can  be  regarded 
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Figure  14:  Probability  distribution  of  the  shock  location  for  the  inhomo¬ 
geneous  Burgers  equations  with  uncertain  initial  condition.  Low  ,  a  =  0.1 
(left)  and  high  input  variability,  a  =  0.6  (right).  Results  obtained  using  the 
Pade-Legendre  stochastic  collocation  approach  [32] 


as  a  preconditioner  (or  a  filter)  such  that  u  ■  Q  is  smooth  and  can  be  efficiently  and 
accurately  interpolated  by  the  polynomial  P. 

Results  obtained  using  the  Pade-Legendre  stochastic  collocation  approach  are 
reported  in  Fig.  5.3  in  terms  of  the  PDF  of  the  shock  locations.  In  this  case  the 
accuracy  is  preserved  even  for  high  input  variability. 

6  Conclusions  and  outlook 

The  growing  interest  in  uncertainty  quantification  has  motivated  the  development 
of  a  variety  of  mathematical  approaches.  In  particular,  probabilistic  uncertainty 
propagation  methods  have  received  considerable  attention. 

The  choice  of  the  appropriate  method  to  use  for  a  specific  application  is  not 
obvious.  For  typical  fluid  mechanics  simulations  some  common  considerations  are: 

•  expensive  function  evaluation :  sampling  based  methods  are  typically  not  ap¬ 
propriate  because  they  might  require  several  thousand  full  computations  to 
build  the  statistics  of  the  outputs; 

•  large  number  of  uncertainties :  boundary  conditions,  material  properties,  ge¬ 
ometry  specification,  etc.  introduce  many  independent  input  parameters  that 
have  to  be  characterize.  Methods  that  suffer  from  curse  of  dimensionality 
quickly  become  unfeasible; 

•  non-linear  system  responses:  transitions  and  bifurcations  are  typical  of  fluid 
mechanics,  especially  for  compressible  flows.  Methods  that  strictly  require  a 
smooth  dependency  between  inputs  and  outputs  can  be  ineffective. 
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In  spite  of  these  difficulties,  probabilistic  UQ  is  naturally  amenable  to  mathe¬ 
matical  and  numerical  analysis  and,  therefore,  is  expected  to  achieve  a  high  level 
of  maturity  in  the  near  future. 

In  addition  to  the  methods  presented  here  several  other  methods  have  been 
applied  especially  in  the  field  of  structural  mechanics.  It  is  also  worth  mentioning 
that  alternative  approaches  not  based  on  probabilistic  reasoning  have  been  proposed 
and  used  with  some  success.  It  is  not  generally  clear  when  probabilistic  methods 
fail  or  are  insufficient;  the  treatment  of  epistemic  uncertainty  remains  difficult  and 
possibly  the  greatest  challenge  in  uncertainty  quantification. 
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